;-------------------------------------------------------------LICENSE--------------------------------------------------------------;
;                                                                                                                                  ;
;The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    ;
;and Message Passing Interface (MPI) parallelization.                                                                              ;
;                                                                                                                                  ;
;Copyright (C) 2012                                                                                                                ;
;Ronglin Jiang                                                                                                                     ;
;rljiang@ssc.net.cn                                                                                                                ;
;585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             ;
;                                                                                                                                  ;
;This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         ;
;as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             ;
;                                                                                                                                  ;
;This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    ;
;MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        ;
;                                                                                                                                  ;
;You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     ;
;Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   ;
;                                                                                                                                  ;
;-------------------------------------------------------------LICENSE--------------------------------------------------------------;

pro map_plot_2d, VAR_NAME = var_name, MESH_FLAG = mesh_flag, CONTOUR_NAME = contour_name, ARROW_NAME = arrow_name, PS = ps,        $
   OUTPUTPATH = outputpath, DATAPATH = datapath, XSIZE = xsize, YSIZE = ysize, RGB_TABLE = rgb_table, VAR_LEVEL = var_level,       $
   CHAR_COLOR = char_color, LINE_COLOR = line_color, CPU_COLOR = cpu_color, BACKGROUND = background, MAG_LV = mag_lv,              $
   PLOT_RANGE = plot_range, ARROW_LEN = arrow_len, ARROW_SCL = arrow_scl, LOG = log, COLORBAR_FORMAT = colorbar_format,            $
   MAG_DIR = mag_dir, MAG_POS_START = mag_pos_start, VAR_OUT = var_out, READDATA_ONLY = readdata_only,                             $
   VAR_MIN = var_min, VAR_MAX = var_max, MAG_PLOT_METHOD = mag_plot_method, OVERPLOT_LV = overplot_lvt, MESH_FILL = mesh_fill,     $
   TICK_LOG = tick_log, COLORBAR_FLAG = colorbar_flag, TIME_START = time_start, TIME_END = time_end, COLOR_TYPE = color_type,      $
   MAG_FIELD_SHIFT = mag_field_shift, TICKS = ticks, ARROW_NX = arrow_nx, ARROW_NY = arrow_ny, SWAP_ENDIAN = swap_endian,          $
   QUIET = quiet

   if keyword_set(VAR_NAME) then begin
      var_flag = 1
   endif else begin
      var_flag = 0
      var_name = 'ro'
   endelse

   if not keyword_set(mesh_flag) then mesh_flag = 0
   if not keyword_set(mesh_fill) then mesh_fill = 0
   if not keyword_set(xsize) then xsize = 1200
   if not keyword_set(ysize) then ysize = 1200
   if not keyword_set(ps) then ps = 0
   if not keyword_set(mag_lv) then mag_lv = 20
   if not keyword_set(log) then log = 0
   if not keyword_set(tick_log) then tick_log = 0
   if not keyword_set(readdata_only) then readdata_only = 0
   if not keyword_set(swap_endian) then swap_endian = 0
   if not keyword_set(quiet) then quiet = 0
   if n_elements(arrow_nx) eq 0 then arrow_nx = 20
   if n_elements(arrow_ny) eq 0 then arrow_ny = 20
   if n_elements(var_level) eq 0 then var_level = 150
   if n_elements(rgb_table) eq 0 then rgb_table = -1
   if n_elements(char_color) eq 0 then char_color = 0
   if n_elements(line_color) eq 0 then line_color = 255
   if n_elements(background) eq 0 then background = 255
   if n_elements(arrow_len) eq 0 then arrow_len = 1.0
   if n_elements(colorbar_format) eq 0 then colorbar_format = '(f4.1)'
   if n_elements(time_start) eq 0 then time_start = 1
   if n_elements(time_end) eq 0 then time_end = 100
   if n_elements(arrow_scl) eq 0 then arrow_scl = 1
   if n_elements(mag_dir) eq 0 then mag_dir = 1
   if n_elements(mag_plot_method) eq 0 then mag_plot_method = 1
   if n_elements(color_type) eq 0 then color_type = 1
   if n_elements(mag_field_shift) eq 0 then mag_field_shift = 0.0
   if n_elements(ticks) eq 0 then ticks = 6

   if quiet eq 0 then begin 
      print
      print, '|-------------------------------------------------------------------------|'
      print, ' Variable name: ', (var_flag eq 1)? var_name : 'none'
      print, ' PS? ', (ps eq 1)? 'Yes' : 'No'
      print, ' Plot mesh? ', (mesh_flag eq 1)? 'Yes' : 'No'
      print, ' Plot contour field: ', (n_elements(contour_name) ne 0)? contour_name : 'NONE'
      print, ' Plot arrow field: ', (n_elements(arrow_name) ne 0)? arrow_name : 'NONE'
      print, '|-------------------------------------------------------------------------|'
      print
   endif

;   var_total_name = [['ro', 'Density'], 
;                     ['gp', 'Gas pressure'], 
;                     ['mp', 'Magnetic pressure'], 
;                     ['te', 'Temperature'], 
;                     ['beta', 'Plasma beta'], 
;                     ['az', 'potential field'], 
;                     ['vx', 'Velocity (x)'], 
;                     ['vy', 'Velocity (y)'], 
;                     ['vz', 'Velocity (z)'], 
;                     ['vv', 'Gas pressure'], 
;                     ['bx', 'Gas pressure'], 
;                     ['by', 'Gas pressure'], 
;                     ['bz', 'Gas pressure'], 
;                     ['bv', 'Gas pressure'], 
;                     ['jx', 'Gas pressure'], 
;                     ['jy', 'Gas pressure'], 
;                     ['jz', 'Gas pressure'], 
;                     ['jv','Gas pressure'], 
;                     ['fgx', 'Gas pressure'], 
;                     ['fgy', 'Gas pressure'], 
;                     ['fgv', 'Gas pressure'], 
;                     ['fggx', 'Gas pressure'], 
;                     ['fggy', 'Gas pressure'], 
;                     ['fggv', 'Gas pressure'], 
;                     ['fmgx', 'Gas pressure'], 
;                     ['fmgy', 'Gas pressure'], 
;                     ['fmgv', 'Gas pressure'], 
;                     ['fmtx', 'Gas pressure'], 
;                     ['fmty', 'Gas pressure'], 
;                     ['fmtz', 'Gas pressure'], 
;                     ['fmtv', 'Gas pressure']]

   usr_var, xsize, ysize
   save_info, entry_info
   usr_font

   if rgb_table eq -1 then begin
      usr_color, color = color_type
   endif else begin
      device, decomposed = 0
      loadct, rgb_table
   endelse

   on_error, 2
   
   if file_test(datapath+'/par.mapvis') eq 0 then begin
      print
      print, ' Error:'
      print, ' Can not find file: par.mapvis'
      print
      retall
   endif

   cd, current = dir_old

;----------------------------------------------------------------------------------------------------------------------------------|
   close, /all

   directory = '' & n_dims = 0 & sys = 0 & ng = 0 & ini_block_x = 0 & ini_block_y = 0 & cpu_number = 0 & output_real8_flag = 0
   nx_no_bnd = 0 & nx_block = 0 & nx_total = 0 & ny_no_bnd = 0 & ny_block = 0 & ny_total = 0 & ratio = 0 & n_levels = 0
   div_clean_flag = 0 & output_level = 0
   nt = 0l & n_var = 0l

   openr, 1, datapath + '/par.mapvis'
   readf, 1, directory
   directory = strcompress(directory, /remove_all)
   readf, 1, cpu_number
   readf, 1, n_dims
   readf, 1, sys
   readf, 1, ng
   readf, 1, ini_block_x
   readf, 1, ini_block_y
   readf, 1, nx_no_bnd
   readf, 1, nx_block
   readf, 1, nx_total
   readf, 1, ny_no_bnd
   readf, 1, ny_block
   readf, 1, ny_total
   readf, 1, ratio
   readf, 1, n_levels
   readf, 1, div_clean_flag
   readf, 1, output_real8_flag
   readf, 1, output_level
   readf, 1, gm
   readf, 1, rangex
   readf, 1, rangey
   readf, 1, zero_posx
   readf, 1, zero_posy
   readf, 1, cfl
   readf, 1, nt
   t = dblarr(nt, /nozero)
   dt = dblarr(nt, /nozero)
   t_tmp = 0.0d0
   for i_t = 0, nt -1  do begin
      readf, 1, t_tmp
      t[i_t] = t_tmp
   endfor
   for i_t = 0, nt -1  do begin
      readf, 1, t_tmp
      dt[i_t] = t_tmp
   endfor
   readf, 1, n_var
   var_names = strarr(n_var)
   var_maxs = fltarr(n_var)
   var_mins = fltarr(n_var)
   var_dims = intarr(n_var)
   var_name_tmp = ''
   for i_var = 0, n_var - 1 do begin
      readf, 1, var_name_tmp
      var_names[i_var] = strcompress(strmid(var_name_tmp, 0, 12), /remove_all)
   endfor
   var_max_tmp = 0.0
   for i_var = 0, n_var - 1 do begin
      readf, 1, var_max_tmp
      var_maxs[i_var] = var_max_tmp
   endfor
   var_min_tmp = 0.0
   for i_var = 0, n_var - 1 do begin
      readf, 1, var_min_tmp
      var_mins[i_var] = var_min_tmp
   endfor
   var_dim_tmp = 0.0
   for i_var = 0, n_var - 1 do begin
      readf, 1, var_dim_tmp
      var_dims[i_var] = var_dim_tmp
      if var_dims[i_var] eq 3 then var_dims[i_var] = 2
   endfor
   close, 1

;----------------------------------------------------------------------------------------------------------------------------------|
   n_blocks = 0l
   dx = 0.0d0
   dy = 0.0d0
   level = 0l
   nc = 0l

   nx_lv = (nx_total - ng * 2) * 2 ^ (output_level - 1) + 1
   ny_lv = (ny_total - ng * 2) * 2 ^ (output_level - 1) + 1

   x = dblarr(nx_no_bnd + 1, /nozero)
   y = dblarr(ny_no_bnd + 1, /nozero)
   x_lv = dblarr(nx_lv, /nozero)
   y_lv = dblarr(ny_lv, /nozero)

   if output_real8_flag eq 8 then begin
      var_tmp = dblarr(nx_no_bnd + 1, ny_no_bnd + 1)
      var_lv_tmp = dblarr(nx_lv, ny_lv)
      vars = dblarr(nx_no_bnd + 1, ny_no_bnd + 1, total(var_dims), /nozero)
      vars_lv = dblarr(nx_lv, ny_lv, total(var_dims), /nozero)
   endif else begin
      var_tmp = fltarr(nx_no_bnd + 1, ny_no_bnd + 1)
      var_lv_tmp = fltarr(nx_lv, ny_lv)
      vars = fltarr(nx_no_bnd + 1, ny_no_bnd + 1, total(var_dims), /nozero)
      vars_lv = fltarr(nx_lv, ny_lv, total(var_dims), /nozero)
   endelse

   if n_elements(outputpath) then cd, outputpath

   if time_start gt nt then begin
      print, ' Warning: The given time_start is too large, which will be the last time.'
      time_start = nt
   endif

   if time_end gt nt then begin
      print, ' Warning: The given time_end is too large, which will be the last time.'
      time_end = nt 
   endif

   if time_start lt 1 then begin
      print, ' Warning: The given time_start is too small, which will be the first time.'
      time_start = 1
   endif
   if time_end lt 1 then begin
      print, ' Warning: The given time_end is too small, which will be the first time.'
      time_end = 1
   endif

;----------------------------------------------------------------------------------------------------------------------------------|
   for ti = time_start, time_end do begin

      datavis_name = datapath + '/data' + name_seq(strcompress(ti, /remove_all), '4.0', /nosymbol) + '.mapvis'
      openr, unit_data, datavis_name, /get_lun, /f77_unformatted, swap_endian = swap_endian

   if quiet eq 0 then print, ' Time = ', t[ti - 1]

      var_id = 0
      if readdata_only eq 0 then begin

         for i_var = 0, n_var - 1 do begin
            if var_name eq var_names[i_var] then begin
               range = [var_mins[i_var], var_maxs[i_var]]
               if n_elements(var_min) ne 0 then range[0] = var_min
               if n_elements(var_max) ne 0 then range[1] = var_max
               title = strupcase(var_name)
               break
            endif
         endfor
         if i_var ge n_var then begin
            message, ' No this variable: ' + var_name
         endif

         var_id = i_var
         if var_flag eq 0 then title = 'novar'
   
         var_levels = findgen(var_level) / (var_level - 1) * (range[1] - range[0]) + range[0]

         xrange = dblarr(2)
         yrange = dblarr(2)
   
         if keyword_set(plot_range) then begin
            xrange[0] = plot_range[0]
            xrange[1] = plot_range[1]
            yrange[0] = plot_range[2]
            yrange[1] = plot_range[3]
         endif else begin
            xrange[0] = 0.0d0 - zero_posx
            xrange[1] = rangex - zero_posx
            yrange[0] = 0.0d0 - zero_posy
            yrange[1] = rangey - zero_posy
         endelse
   
         file = strcompress(title, /remove_all) + '_' +                                                                            $
               'contour' + strcompress(keyword_set(contour_name), /remove_all) + '_' +                                             $
               'arrow' + strcompress(keyword_set(arrow_name), /remove_all) + '_' +                                                 $
               'mesh' + strcompress(mesh_flag, /remove_all) + '_' +                                                                $
               'zoom' + strcompress(keyword_set(plot_range), /remove_all) + '_' +                                                  $
               'time' + name_seq(strcompress(t[ti - 1], /remove_all), '3.3')
         usr_plot_sys, entry_dev, image_name = file, xsize = xsize, ysize = ysize, /start, /pixmap, ps = ps

         pos_x_shift = 0.0
         if not keyword_set(colorbar_flag) then pos_x_shift = 0.05
         contour, var_lv_tmp, x_lv, y_lv, xstyle = 1, ystyle = 1, levels = var_levels,                                             $
            posi = [0.15 + pos_x_shift, 0.20, 0.80 + pos_x_shift, 0.85],                                                           $
            background = background, color = char_color, /fill, title = 'Time=' + strmid(strcompress(t[ti - 1]), 1, 4),                $
            xtitle = 'x', ytitle = 'y', nodata = 1, xrange = xrange, yrange = yrange,                                              $
            charsize = !usr_chars, thick = !usr_thick, charthick = !usr_chart, xthick = !usr_xthick, ythick = !usr_ythick
   
         if keyword_set(colorbar_flag) and tick_log eq 1 then colorbar, posi = [0.82, 0.20, 0.86, 0.60], /vertical,                $
            title = title + ' (log)', /right, range = alog10(range), format = colorbar_format, background = background,            $
            color = char_color, charsize = !usr_chars_bar, charthick = !usr_chart
         if keyword_set(colorbar_flag) and tick_log eq 0 then colorbar, posi = [0.82, 0.20, 0.86, 0.60], /vertical,                $
            title = title, /right, range = range, format = colorbar_format, background = background, color = char_color,           $ 
            ylog = log, yticks = ticks, charsize = !usr_chars_bar, charthick = !usr_chart
      endif

;----------------------------------------------------------------------------------------------------------------------------------|
      for cpui = 0, cpu_number - 1 do begin

         readu, unit_data, n_blocks

         for i_block = 0, n_blocks - 1 do begin

            readu, unit_data, level, nc
            readu, unit_data, dx, dy
            readu, unit_data, x
            readu, unit_data, y

            for i_var = 0, total(var_dims) - 1 do begin
               readu, unit_data, var_tmp
               vars[*, *, i_var] = var_tmp
            endfor

            if nc eq 0 then begin
               if var_dims[var_id] eq 2 then begin
                  var = sqrt(vars[*, *, total(var_dims[0:var_id])-1] ^ 2 + vars[*, *, total(var_dims[0:var_id])-2] ^ 2)
               endif else begin
                  var = vars[*, *, total(var_dims[0:var_id])-1]
               endelse

               if (keyword_set(var_max)) then begin
                  tmp_index = where(var gt var_max)
                  if n_elements(tmp_index) ne 1 then var[tmp_index] = var_max
               endif
               if (keyword_set(var_min)) then begin
                  tmp_index = where(var lt var_min)
                  if n_elements(tmp_index) ne 1 then var[tmp_index] = var_min
               endif

               if keyword_set(log) then begin
                  var = alog10(var)
                  range_log = alog10(range)
                  var_levels = findgen(var_level) / (var_level - 1) * (range_log[1] - range_log[0]) + range_log[0]
               endif

;----------------------------------------------------------------------------------------------------------------------------------|
               if readdata_only eq 0 then contour, var, x, y, levels = var_levels, /overplot, /fill, nodata = 1 - var_flag

            endif

            if readdata_only eq 0 then begin
               xb1 = min(x)
               xb2 = max(x)
               yb1 = min(y)
               yb2 = max(y)
   
               if level eq -1 then polyfill, [xb1, xb2, xb2, xb1, xb1], [yb1, yb1, yb2, yb2, yb1], color = char_color
   
               if keyword_set(cpu_color) then begin
                  if mesh_fill eq 1 then polyfill, [xb1, xb2, xb2, xb1, xb1], [yb1, yb1, yb2, yb2, yb1],                    $ 
                     color = 240 * cpui / cpu_number
;                     color =  240 * (cpui / 4) / (cpu_number / 4)
                  if mesh_flag eq 1 then oplot, [xb1, xb2, xb2, xb1, xb1], [yb1, yb1, yb2, yb2, yb1],                       $
                     thick = round(!usr_thick / 1.5), color =  240 * cpui / cpu_number
               endif else begin
                  if mesh_fill eq 1 then polyfill, [xb1, xb2, xb2, xb1, xb1], [yb1, yb1, yb2, yb2, yb1],                    $ 
                     color = 240 * cpui / cpu_number
                  if mesh_flag eq 1 then oplot, [xb1, xb2, xb2, xb1, xb1], [yb1, yb1, yb2, yb2, yb1],                       $
                     thick = round(!usr_thick / 1.5), color = line_color
               endelse
            endif
         endfor
      endfor

      if readdata_only eq 0 then begin
         axis, xaxis = 0, xstyle = 1, color = char_color, xtick_get = x_tick,                                                      $
            charsize = !usr_chars, charthick = !usr_chart, xthick = !usr_xthick, ythick = !usr_ythick
         axis, yaxis = 0, ystyle = 1, color = char_color, ytick_get = y_tick,                                                      $
            charsize = !usr_chars, charthick = !usr_chart, xthick = !usr_xthick, ythick = !usr_ythick
         x_tick_len = n_elements(x_tick)
         y_tick_len = n_elements(y_tick)
         axis, xaxis = 1, xstyle = 1, color = char_color, xtickname = replicate(' ', x_tick_len),                                  $
            charsize = !usr_chars, charthick = !usr_chart, xthick = !usr_xthick, ythick = !usr_ythick
         axis, yaxis = 1, ystyle = 1, color = char_color, ytickname = replicate(' ', y_tick_len),                                  $
            charsize = !usr_chars, charthick = !usr_chart, xthick = !usr_xthick, ythick = !usr_ythick
      endif

      readu, unit_data, level, nc
      readu, unit_data, dx, dy
      readu, unit_data, x_lv
      readu, unit_data, y_lv

      for i_var = 0, total(var_dims) - 1 do begin
         readu, unit_data, var_lv_tmp
         vars_lv[*, *, i_var] = var_lv_tmp
      endfor

;----------------------------------------------------------------------------------------------------------------------------------|
;     var output
;----------------------------------------------------------------------------------------------------------------------------------|
      var_out = create_struct('header', 'This is created by map_plot_2d of MAP code.')
      var_out = create_struct(var_out, 't', t[ti - 1])
      var_out = create_struct(var_out, 'x', x_lv, 'y', y_lv)
      var_out = create_struct(var_out, 'nx', nx_lv, 'ny', ny_lv)
      for i_var = 0, n_var - 1 do begin
         if var_dims[i_var] eq 2 then begin
            var_out = create_struct(var_out, var_names[i_var], vars_lv[*, *, floor(total(var_dims[0:i_var])-2) :              $
                                                                             floor(total(var_dims[0:i_var])-1)])
         endif else begin
            var_out = create_struct(var_out, var_names[i_var], vars_lv[*, *, floor(total(var_dims[0:i_var])-1)])
         endelse
      endfor

      if readdata_only ne 0 then continue

;----------------------------------------------------------------------------------------------------------------------------------|
;     overplot contour
;----------------------------------------------------------------------------------------------------------------------------------|
      if n_elements(overplot_lv) ne 0 then begin
         var = vars_lv[*, *, var_id]

         annotation = replicate ('', n_elements(overplot_lv))
         for i = 0, n_elements(overplot_lv) - 1 do annotation[i] = name_seq (overplot_lv[i], '3.2')
         if keyword_set(log) then begin
            var = alog10(var)
            contour, var, x_lv, y_lv, levels = alog10(overplot_lv), /overplot, nodata = 1 - var_flag,                              $
               thick = round(!usr_thick / 2), c_labels = replicate(1, n_elements(overplot_lv)),                                    $
               c_charsize = !usr_chars / 2., c_charthick = !usr_chart, c_annotation = annotation
         endif else begin
            contour, var, x_lv, y_lv, levels = overplot_lv, /overplot, nodata = 1 - var_flag,                                      $
               thick = round(!usr_thick / 2), c_labels = replicate(1, n_elements(overplot_lv)),                                    $
               c_charsize = !usr_chars / 2., c_charthick = !usr_chart, c_annotation = annotation
         endelse
      endif

;----------------------------------------------------------------------------------------------------------------------------------|
;     contour field
;----------------------------------------------------------------------------------------------------------------------------------|
      if keyword_set(contour_name) eq 1 then begin
         for i_var = 0, n_var - 1 do begin
            if contour_name eq var_names[i_var] and var_dims[i_var] eq 1 then break
         endfor
         if i_var ge n_var then begin
            message, ' No this variable: ', contour_name
         endif

         eps = 0.0
         if var_maxs[i_var] le var_mins[i_var] then eps = 1e-6
         if mag_plot_method eq 1 then begin
            levels = findgen(mag_lv + 1) / mag_lv * (var_maxs[i_var] - var_mins[i_var] + eps) + var_mins[i_var] +                  $
               mag_field_shift * (var_maxs[i_var] - var_mins[i_var] + eps) / mag_lv

            contour, vars_lv[*, *, total(var_dims[0:i_var]) - 1], x_lv, y_lv,                                                      $
               color = line_color, /overplot, levels = levels, thick = !usr_othick * 1.5
         endif else begin
;            vector_potential_2d, bx_lv, by_lv, x_lv, y_lv, kesai, dir = mag_dir, pos_start = mag_pos_start

;            max_k = max(kesai) + eps
;            min_k = min(kesai)
;            levels = findgen(mag_lv + 1) / mag_lv * (max_k - min_k) + min_k

;            contour, kesai, x_lv, y_lv, color = line_color, /overplot, levels = levels, thick = !usr_othick
         endelse
      endif

;----------------------------------------------------------------------------------------------------------------------------------|
;     arrow
;----------------------------------------------------------------------------------------------------------------------------------|
      if keyword_set(arrow_name) eq 1 then begin
         for i_var = 0, n_var - 1 do begin
            if arrow_name eq var_names[i_var] and var_dims[i_var] eq 2 then break
         endfor
         if i_var ge n_var then begin
            message, ' No this variable: ', contour_name
         endif

         vx_lv = vars_lv[*, *, total(var_dims[0:i_var]) - 2]
         vy_lv = vars_lv[*, *, total(var_dims[0:i_var]) - 1]

         c = 10 / rangey * (xrange[1] - xrange[0] + yrange[1] - yrange[0]) / 2 / arrow_len *                                       $
            (rangex + rangey) / (xrange[1] - xrange[0] + yrange[1] - yrange[0])

         devmax = convert_coord([0, 0], [0, (xrange[1] - xrange[0] + yrange[1] - yrange[0]) / c * 1.0], /data, /to_dev)
         lenmax = sqrt((devmax[0, 1] - devmax[0, 0]) ^ 2 + (devmax[1, 1] - devmax[1, 0]) ^ 2) / 2.5
         for ii = 0, arrow_nx - 1 do begin
            a0 = (ii + 0.5) * (xrange[1] - xrange[0]) / arrow_nx + xrange[0]
            for jj = 0, arrow_ny - 1 do begin
               b0 = (jj + 0.5) * (yrange[1] - yrange[0]) / arrow_ny + yrange[0]

               for i = 0, nx_lv - 1 do if x_lv[i] ge a0 then break
               for j = 0, ny_lv - 1 do if y_lv[j] ge b0 then break

               vx_arrow = vx_lv[i, j] * (a0 - x_lv[i - 1]) * (b0 - y_lv[j - 1]) +                                                  $
                          vx_lv[i - 1, j] * (x_lv[i] - a0) * (b0 - y_lv[j - 1]) +                                                  $
                          vx_lv[i - 1, j - 1] * (x_lv[i] - a0) * (y_lv[j] - b0) +                                                  $
                          vx_lv[i, j - 1] * (a0 - x_lv[i - 1]) * (y_lv[j] - b0)
               vy_arrow = vy_lv[i, j] * (a0 - x_lv[i - 1]) * (b0 - y_lv[j - 1]) +                                                  $
                          vy_lv[i - 1, j] * (x_lv[i] - a0) * (b0 - y_lv[j - 1]) +                                                  $
                          vy_lv[i - 1, j - 1] * (x_lv[i] - a0) * (y_lv[j] - b0) +                                                  $
                          vy_lv[i, j - 1] * (a0 - x_lv[i - 1]) * (y_lv[j] - b0)

               vx_arrow = vx_arrow / (x_lv[i] - x_lv[i - 1]) / (y_lv[j] - y_lv[j - 1])
               vy_arrow = vy_arrow / (x_lv[i] - x_lv[i - 1]) / (y_lv[j] - y_lv[j - 1])

               a1 = a0 + (xrange[1] - xrange[0] + yrange[1] - yrange[0]) / c * vx_arrow / arrow_scl
               b1 = b0 + (xrange[1] - xrange[0] + yrange[1] - yrange[0]) / c * vy_arrow / arrow_scl
               devp = convert_coord([a0, a1], [b0, b1], /data, /to_dev)
               hsize = sqrt((devp[0, 1] - devp[0, 0]) ^ 2 + (devp[1, 1] - devp[1, 0]) ^ 2) / 2.0
               if hsize gt lenmax then hsize = lenmax

               arrow, a0, b0, a1, b1, color = line_color, hsize = hsize, /data, /solid, thick = !usr_othick * 1.6
            endfor
         endfor

         a1 = xrange[1] + (xrange[1] - xrange[0]) / 12.0
         b1 = yrange[1] - (yrange[1] - yrange[0]) / 5.0 + (xrange[1] - xrange[0] + yrange[1] - yrange[0]) / c * 1.0
         hsize = lenmax
         arrow, xrange[1] + (xrange[1] - xrange[0]) / 12.0, yrange[1] - (yrange[1] - yrange[0]) / 5.0,                             $
            a1, b1, color = 255 - line_color, hsize = hsize, /data, /SOLID, thick = !usr_othick
         xyouts, xrange[1] + (xrange[1] - xrange[0]) / 15.0, yrange[1] - (yrange[1] - yrange[0]) / 4.0,                            $
            '!4v!3 = ' + strcompress(arrow_scl, /remove_all), color = 255 - line_color,                                            $
            charsize = !usr_chars_bar , charthick = !usr_chart
      endif

;----------------------------------------------------------------------------------------------------------------------------------|
      usr_plot_sys, entry_dev, image_name = file, /done, ps = ps
      close, unit_data
      free_lun, unit_data
   endfor

;----------------------------------------------------------------------------------------------------------------------------------|
   cd, dir_old
   reset_info, entry_info

end

